Computational study of subcritical response in flow past a circular cylinder 
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Flow past a circular cylinder is investigated in the subcritical regime, below the onset of Benard- 
von Karman vortex shedding at Rec — 47. The transient response of infinitesimal perturbations is 
computed. The domain requirements for obtaining converged results is discussed at length. It is 
shown that energy amplification occurs as low as Re = 2.2. Throughout much of the subcritical 
regime the maximum energy amplification increases approximately exponentially in the square of Re 
reaching 6800 at Rec. The spatiotemporal structure of the optimal transient dynamics is shown to 
be transitory Benard-von Karman vortex streets. At Re ~ 42 the long-time structure switches from 
exponentially increasing downstream to exponentially decaying downstream. Three-dimensional 
computations show that two-dimensional structures dominate the energy growth except at short 
times. 

PACS numbers: 47.20.Ft, 47.15. Tr, 47.10. ad, 47.11. Kb 
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I. INTRODUCTION 

Incompressible fluid flow past a circular cylinder has 
been extensively studied, both for its relevance to nu- 
merous engineering applications and as a prototype bluff- 
body flow exhibiting vortex shedding. See for exam- 
ple [HH] and references therein. It is one of the most-used 
test bed for exploring stability concepts in open flows, 
e.g. f3HT3). As a result, a great deal is known about this 
flow in general and in particular concerning the primary 
instability. It is well-established that this instability oc- 
curs at a critical Reynolds number of about 47 [SHl]. 
Below this value the steady wake flow is linearly stable, 
while above it the steady flow is unstable and periodic os- 
cillations arise leading to the famous Benard-von Karman 
vortex street [Ml [15] . Our concern here is what happens 
in the stable, subcritical regime prior to, and leading up 
to, the onset of oscillations. 

Stable flows may exhibit transient growth (THl ITT] . 
This means that inflnitesimal perturbations to the flow 
may grow in energy for some time before subsequently de- 
caying to zero. While initially popular in parallel shear 
flows as possibly playing a role in the transition to turbu- 
lence, e.g. [I8H20] . transient growth has become increas- 
ingly of interest in spatially developing flows, e.g. pTll^TI - 
125] . For instance, separated flows arising due to abrupt 
changes in geometry are known to promote extremely 
large transient growth in perturbations ^31 [26l [27] . The 
origin of this growth can be traced to the non-normality 
of the linear stability operator associated with many 
shear flows [TT|[16]. This means, in particular, that in 
spatially developing flows the eigenmodes of the stabil- 
ity operator tend to be located downstream while the 
eigenmodes of the adjoint operator tend to be located 
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upstream [TTl [15] [TT] . 

For the cylinder wake, Giannetti and Luchini [T3] first 
examined in detail the adjoint eigenmodes in the vicin- 
ity of the primary instability and used these, together 
with direct eigenmodes, to understand the sensitivity of 
the flow. Their results are further discussed in detail 
by Chomaz in the context of non-normality [TT]. Since 
this important work, there have been further computa- 
tions of direct and adjoint modes and transient growth 
for the cylinder wake. For example Marquet et al. [25] 
have computed direct and adjoint eigenmodes of the su- 
percritical flow, and Abdessemed et al. [29J have studied 
the transient growth, focusing on supercritical Reynolds 
numbers, although also reporting some subcritical values. 

There have been a number of experimental studies of 
the cylinder wake in the stable and marginally unstable 
regime [3] [301 [SI]- The most relevant are the studies 
by Le Gal and Croquette [30] and the recent work by 
Marais et al. [3T] on the impulse response at subcritical 
Reynolds numbers. By inducing an impulse through a 
small displacement or rotation of the cylinder, wavepack- 
ets are generated that grow and are subsequently ad- 
vected downstream. While the measurements by Le Gal 
and Croquette provide informative qualitative proper- 
ties of the transient dynamics, these measurements were 
made using streaklines and so provide limited quanti- 
tative detail. The more recent work by Marais et al. 
uses particle image velocimetry to obtain quantitative 
measurements of the transient response in the subcriti- 
cal regime. 

The purpose of the current study is two-fold. Primar- 
ily we establish an accurate characterization of the opti- 
mal transient energy growth throughout the subcritical 
regime for the cylinder wake. We determine the thresh- 
old Reynolds number where energy growth first occurs, 
determine the Reynolds number dependence of the opti- 
mal growth, and its value at criticality We show that 
the transient dynamics associated with optimal energy 
growth is in the form of wave packets similar to those 
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FIG. 1. Diagram of the cylinder geometry (not to scale), 
showing the inflow, outflow and cross-stream dimensions ref- 
erenced later. Also marked are the separation streamlines and 
the downstream stagnation point, a;^. 



observed in experiments on subcritical wakes. The sec- 
ondary purpose of the paper is to highlight and establish 
the computational requirements for such computations. 
As we shall show, the requirements for accurate compu- 
tations of transient growth are more severe than those of 
linear stability. While these findings are specific to the 
cylinder wake, they should guide computations of similar 
flows. 



II. FORMULATION 

The flow geometry is illustrated in Fig. [T] A circular 
cylinder of diameter D is placed in a free-stream flow 
Uco- Streamwise x and cross-stream y coordinates are 
centered on the circular cross section. The cylinder axis, 
infinite in length and normal to the free-stream velocity, 
aligns with the z-coordinate. 

In principle this open flow would have infinite extent 
in all directions. In practice, however, our numerical cal- 
culations necessarily employ a computational domain f2 
with finite inflow Li, outflow io, and cross-stream Lc 
lengths, as illustrated. The z-direction is homogeneous, 
and for the issues addressed in this paper, this direction 
can be treated without needing to restrict to a bounded 
domain. The demands on the domain dimensions is 
an important aspect of our work discussed in detail in 

secnfn 

The fluid is governed by the incompressible Navier- 
Stokes equations 



dtn+ (u- V)u = -Vp- 
V • u = 0, 



Re-^V^u, 



(la) 
(lb) 



where u — u(x) = {u{x, y, z), v{x, y, z), w{x, y, z)) is the 
fluid velocity and p(x,y,z) is the static pressure. With- 
out loss of generality we set the density to unity. The 



equations are non-dimensionalized by the free-stream 
speed Uoo and the cylinder diameter D. The Reynolds 
number Re is therefore given as 



Re: 



U^D 



where v is the kinematic viscosity of the fluid. 

No-slip boundary conditions are imposed on the cylin- 
der surface. The boundary conditions around the outer 
boundaries of the domain are such as to give a good nu- 
merical approximation of the unbounded flow. Specifl- 
cally, the boundary conditions are: 

u{dQ.^,t) = Uooex, (2a) 

U(90c,i) = ^ooCa;, (2b) 

u{dn^,t)^{), (2c) 

dMd^o, t) - 0, p{d^o. t) = 0, (2d) 

where d^i is the inlet boundary at a; = — ii, 9ilc is the 
cross-stream boundary at y = ±Lc, dflu, is the boundary 
of the cylinder, and d^o is the outlet boundary at x = Lo- 

The remaining material in this section is included for 
completeness and to clearly define notation. Since the 
details are contained in numerous prior publications, es- 
pecially [32l[33], the treatment here is minimal. 

Equation M is solved using Direct Numerical Simu- 
lation (DNS) employing a split-step pressure-correction 
scheme described elsewhere [Ml |3S] . This is implemented 
in a spectral-element code [36] utilizing an elemental de- 
composition of the domain in the two-dimensional (2D) 
plane normal to the cylinder axis. 

The base flows U considered in this paper are steady, 
two-dimensional solutions to Eq. (fTl). Hence U — 
{U{x,y),V{x,y)). These Reynolds number dependent 
flows are symmetric about the streamwise centerline as 
depicted in Fig. [l] Figure [2| shows a typical base flow. 
Those at other Re are qualitatively similar, differing pri- 
marily in the length of the recirculation region behind the 
cylinder. For Re < 6.2, there is no recirculation region. 
Steady base flows in both subcritical and supercritical 
regimes are rapidly obtained through DNS by imposing 
this midplane symmetry. Once computed, base flows are 
stored for use in subsequent linear calculations. 

Our interest is in the dynamics of infinitesimal pertur- 
bations u' to the steady base flow. These perturbations 
evolve according to the linearized Navier-Stokes equa- 
tions 



dtu' + (u' • V)U + (U • V)u' = -Vp' 



0, 



Rc-iyV 



(3a) 
(3b) 



where p' is the perturbation pressure. Numerically 
Eq. ([3]) is solved using the same techniques as the non- 
linear Navier-Stokes equations. For the most part we 
will focus on 2D perturbation fields u' = {u',v') on 
2D grids. However, we will also consider briefly three- 
dimensional perturbations. Since the base flow is 2D, 




FIG. 2. Representative base flow. Streamlines of the flow at 
Re=40. The computational domain is much larger than the 
region shown. The main flow field uses a contour spacing of 
0.25, while a smaller spacing of 0.01 is used to highlight the 
structure of the recirculation bubble. 



three-dimensional (3D) perturbations can be decomposed 
into non-interacting modes of the form 



u'{x, y, z) — u{x, 2/)e*^^ -I- a 



(4) 



where /3 is the spanwise wavenumber of the perturba- 
tion. Only u(x,y), a three-component field on a two- 
dimensional grid, needs to be computed. 

Our primary focus is on the transient dynamics of per- 
turbations at subcritical Reynolds numbers. We focus on 
the energy of perturbation fields and seek initial condi- 
tions u'(0) which generate the largest possible growth in 
energy under evolution by Eq. (|3|. The formalism is as 
follows. Let A{t) denote the linear evolution operator 
over a time t defined by Eq. (Is]), so that 

u\t)=A{t)u'{0). 

Since the governing equations are linear it is sufficient 
to consider initial perturbation fields with ||u'(0)|p — 
(u'(0),u'(0)) = 1, where (•, •) denotes the L2 inner- 
product. Then the energy growth in the perturbation 
field over a time r is given by 



E{i 



EiO) 



= (u'(r),u'(r)). 



In terms of the operator A(t), and its adjoint A*{t) in 
the L2 inner-product, we have 



E{r) 
E{0) 



--{A{r)n'{0),A{T)u'm 
= (u'(0),X(r)^(r)u'(0)). 



Letting Xj and Vj denote eigenvalues and normalized 
eigenfunctions of the operator A*{t)A{t), we have 



A*{T)A{T)^^J 



h^3^ 



1. 



(5) 



The eigenvalues are non-negative and we assume ordering 
Ai > A2 > •••. 



The maximum possible energy growth, denoted G{t), 
over a specified time horizon t, is then given by the dom- 
inant eigenvalue of ^*(t)^(t), i.e. 



G(r) 



max Xj — Ai , 

j 



The initial perturbation leading to this growth is the cor- 
responding eigenfunction Vi. While the dominant eigen- 
value of A*{t)A{t) is generally of most importance, the 
first few sub-dominant eigenvalues may also be of inter- 
est. In particular V2 will also be considered in this study. 
The maximum energy growth over all time horizons is 
denoted by 



G" 



where 



G(r'"'^'^), 



arg max G'(t). 



(6) 



(7) 



While our primary focus is transient growth, we report 
some eigenvalue results. Equation ra can formally be 
written 

dtu' = £u'. 

Looking for normal-mode solutions to these equations 
gives the eigenvalue problem 



Cuj = a 



3"-3' 



\\u,\\ = l, (8) 

where Uj are normalized cigcnmodes and aj eigenvalues 
of C. We assume ordering such that Re(cri) > Re(cr2) > 
• • • . Stability of the base flow is determined from the 
right-most eigenvalues of C in the complex plane. 

Associated to Eq. (l8| is the adjoint eigenvalue problem 



C*u* 



(9) 



where u* are the adjoint modes, (eigenmodes of the ad- 
joint operator £*), and a* is the complex conjugate of 
aj . The norm of the adjoint eigenmodes is chosen so that 
{u*j,Uj) = 1 for all j. Then the eigenmode and adjoint 
modes satisfy biorthonormality: 

{u*,Uj) =S,j. 

In practice the eigenvalues for the transient growth 
or stability problems are computed through a modified 
Arnoldi algorithm using a time-stepper approach [321133) . 



III. INFLUENCE OF DOMAIN SIZE 

As noted in the introduction, the size of the computa- 
tional domain can be an important factor in studies of 
the transient response in subcritical cylinder fiow. While 
the requirements for accurate base flows and eigenvalue 
calculations for the cylinder wake have been discussed in 
many places _^37-40 , and are presented in our study in 
the Appendix, there is no such discussion for transient 
growth calculations for the cylinder wake. Hence some 
details are worthwhile. We first present results from the 
convergence study and then discuss some of the causes 
and implications of our findings. 
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FIG. 3. Representative spectral-elemental mesh. Dimensions 
are Li — 45, Lc = 45 and Lo — 125 (refer to Fig. [l] for 
definitions). 



A. Convergence 

We focus on the role of inflow length, Li, and the 
cross-stream half-length, Lc, since these are the critical 
lengths. The requirements on the outflow length, Lo, 
are set by the largest r value under consideration in the 
transient growth analysis and could in principle be arbi- 
trarily large. Based on the maximum value of r = 110 we 
consider, and a free-stream t/oo ~ 1, we fix the outflow 
length in the convergence study at Lq — 125. 

Figure |3] shows a representative spectral-element do- 
main of the type used in our computations. (It is in 
fact the final mesh used for obtaining results presented 
in Sec. |IV[ ) For the study of domain size, Li and Lc are 
varied by adding or removing elements as necessary, and 
the polynomial order of the spectral expansion within 
each element is fixed at order 8. The polynomial order 
used in obtaining the final transient growth results is 6, 
as established in the Appendix. 

The dependence of these calculations on domain size 
is assessed through the calculation of energy growth at a 
fixed time horizon of r = 20. At this time horizon non- 
negligible growth is expected across most of the range of 
Reynolds numbers under consideration. Figure l4] sum- 
marizes the errors introduced through domain size re- 
striction. Three Reynolds numbers. Re = 5, Re = 20 
and Re = 46, are considered to ensure that the mesh is 
capable of resolving all solutions in the subcritical range. 

The transient growth results are seen to be sensitive 
to domain size, much more so than either the base flow 
or eigenvalue calculations presented in the appendix. Do- 
mains that provide results accurate to within 1% for base 
flows and eigenvalues, e.g. a domain with Li = Lc = 25, 
do not provide such accuracy for transient growth calcu- 
lations. The effect of cross-stream restriction at low Re 
is particularly significant. Even accepting that in many 
cases one does not need high accuracy in transient growth 
values. Fig. |4] demonstrates the care that must be taken 
in computing transient response in the subcritical regime. 

Based on these results, a computational domain with 




FIG. 4. Convergence of optimal growth calculation with mesh 
geometry by (a) in-flow length Li and (b) cross-stream length 
Lc- Points indicate the computed values. Optimal growth is 
for a time horizon of r = 20. Percentage errors are shown 
relative to the calculation using Li — 65 and Lc = 65, respec- 
tively. 



Li ~ 45 and Lc = 45 is deemed sufficient to resolve tran- 
sient growth calculations to within about 1% for subcriti- 
cal Reynolds numbers. Possibly the accuracy is not quite 
1% at Re < 5, but the growth values are sufficient for our 
purposes. A diagram of the resulting mesh is shown in 
Fig. HI 

B. Discussion 

We begin by recalling that recently, Abdessemed et 
al. |29| reported transient growth calculations in the 
cylinder wake, including some within the subcritical 
regime at Re = 45. In comparing those results with 
ours, we have found that their growth values are about 
32% larger than those computed on the mesh in Fig. [3J 
at the same Re and time horizon. We will use this dis- 
crepancy to focus the present discussion. 

The Abdessemed et al. calculations were performed 
using a spectral-element code similar to that used in 
this study. Their computational domain has bounds 
—8 < a; < 95 and —12.5 < y < 12.5 and identical bound- 
ary conditions to ours. One can quickly rule out the 
possibility that resolution (polynomial order) or outflow 
length are significant factors in the disagreement between 
the two computations. Moreover, we have already seen 
that the transient growth calculations require large inflow 
Li and cross-stream Lc dimensions so the disagreement 



is not surprising in retrospect. However, there are in fact 
two causes for the discrepancy which we address: one is 
the indirect effect caused by differences in the base flows 
for the different computations and the other is the di- 
rect effect of domain requirements for the optimal initial 
condition itself. 

We shall refer to our computational domain with di- 
mensions as in Fig. p] as ilj^, (large domain), and that 
with dimensions used by Abdessemed et al. as ilg i (small 
domain). Let Ul denote the base flow computed on IIl 
(at Re = 45) and let u'l denote the normalized initial 
condition, ||u'l|| — 1, giving optimal growth at t = 100. 
Similarly let Us and u's denote the base flow and normal- 
ized optimal initial condition on D.a, . The resulting energy 
growth for the two calculations is given in the first two 
rows of Table |lj where one sees the large discrepancy. It 
is worth pointing out that the critical Reynolds numbers 
obtained on the two meshes differ by only about 2%. 

To assess the role of the base flow, one can take the 
initial condition u'g from the smaller domain and evolve 
it on the larger domain with the corresponding base flow 
Ul. The resulting growth after 100 time units is given in 
the third row of Table |lj Necessarily the growth had to 
be less than for u'l because u'l, by deflnition, gives the 
largest possible growth over this time horizon on JIl • It is 
perhaps somewhat surprising, however, that the growth 
following from the fixed initial condition u'g is approxi- 
mately factor of 2 less on the large domain than on the 
small domain (second and third rows of Table IT]) . The 
difference is almost entirely attributable to the difference 
in the base flows Ul and Ug. This is confirmed by evolv- 
ing u'g on the small domain but with the base flow Ul, 
truncated onto the smaller domain. The result is given in 
the last row of Table HI There is little difference between 
the evolution of u'g on the two domains, if they both 
have the same base flow Ul- The conclusion is that the 
energy growth may depend considerably on the base flow 
(a factor of 2 in this case), even in situations where other 
measures, such as critical Reynolds numbers, would not 
reveal such a large dependency. 

There is then the remaining issue of how u'l and u'g 
differ and why the energy growth following from u'g is 
27% less than from u'l for the same base flow (flrst and 
third rows in Table IT]). This has to do with the domain 
requirements, in particular the inflow length Li needed to 
capture the optimal initial condition. In Fig. [5] we show 
the upstream extent of u'l on two scales. The energy of 
the perturbation upstream of x = —5 is of the order 10^^ 
and, while one might consider it to be negligible, this por- 
tion of the initial condition makes a signiflcant contribu- 
tion to the overall growth and cannot be neglected in the 
transient growth computations if quantitative accuracy 
is required. 

We conclude with a few further remarks on the pres- 
ence of weak upstream tails in the optimal initial condi- 
tions. First, despite our caution about the need to re- 
solve these to obtain quantitatively accurate results, we 
flnd the linear evolution from the optimal initial condi- 
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FIG. 5. (color online) Energy of the optimal initial condition 
for Re = 45, r = 100 shown on two scales: (a) matching 
that of the figures shown later in the paper and, (b) showing 
an extended scale which highlights the upstream tail of the 
perturbation. 



domain 


base flow 


IC 


£(100) 


f^L 


Ul 


u'l 


2.54 X 10^ 


fis 


Us 


u's 


3.36 X 10^ 


t^L 


Ul 


u's 


1.85 X 10^ 


fis 


Ul 


u's 


1.77 X 10^ 



TABLE I. Effects of domain size, base flow, and initial con- 
dition (IC) on the energy growth at Re = 45. Energy at time 
100 is given for different configurations (see text). 



tion is qualitatively similar whether or not the numerical 
domain fully contains the weak upstream tail of the initial 
condition. We observe no important qualitative errors in 
discounting it, but quantitatively the errors in the en- 
ergy growth can be large. In addition, the length of the 
upstream tail depends on the time horizon r. The value 
T=100 considered in our comparison is rather large. For 
smaller time horizons the weak tail may be absent from 
the optimal initial condition simply because such a tail 
could not advect downstream and come into play over a 
small time horizon. Specifically, in Sec. IVB| we focus on 
optimal initial conditions computed for r=20 and in this 
case the upstream tails are absent. 



IV. TRANSIENT SUBCRITIAL RESPONSE 

A. 2D Energy Growth 

Figures |6] and [7] summarize the optimal energy growth 
for 2D perturbations in the subcritical regime. Figure [6] 
shows the optimal growth envelopes for particular values 
of Re. To be clear, these curves show the largest attain- 
able energy growth over all possible initial conditions at 
each value of r. The uppermost curve is the growth at 
Re = 50, above the onset of linear instability at Rcc. Af- 
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FIG. 6. Optimal energy growth at Reynolds numbers from 
Re = 5 to Re = 50 in increments of Re = 5. Points indicate 
the computed values. The case Re = 50 is above Rec and is 
shown as dashed. The horizontal line is an estimate of the 
maximum growth in the subcritical regime (see text). 



FIG. 8. Maximum growth C""^ as a function of the square 
of Reynolds number. The vertical dotted line corresponds 
to Rcc. G™"^, the maximum growth at Rec, is indicated 
with a point and is computed differently from cases Re < 46 
(see text). The dashed-dotted line highlights the relationship 
given by Eq. (10 1. 
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FIG. 7. Contour plot of optimal energy growth in the subcrit- 
ical regime, with contour levels as indicated. The thick black 
curve denotes the contour of no-growth, G — 1. 



ter an initial rapid growth, the energy increase saturates 
to an exponential rate, in line with that of the leading 
eigenvalue. 

Figure [t] shows growth contours in the (Re, r) plane. 
The contours highlight the fast energy growth at small 
time horizons and the slow decay for long r. The thick 
curve denotes the no-growth contour: G = 1. The inter- 
ception of this curve with the Re-axis indicates a criti- 
cal Reynolds number for energy growth [41 , Reg, below 
which all perturbations decay monotonically in time and 
above which there is transient energy growth for at least 
some perturbations. We estimate Re^ = 2.2. Thus, a 
small amount of transient energy growth is possible be- 
fore the formation of the recirculation region behind the 
cylinder at Re = 6.2. 



The single most important measure of the transient en- 
ergy growth at any Re is the maximum G"^^^ over all r. 
[Recall Eqs. ^ and ^]. This is shown in Fig. [s] where 
logjg G'^^^ is plotted as a function of Re^. Throughout 
most of the subcritical regime, the maximum growth in- 
creases exponentially with Re . More specifically, we find 



G" 



exp(1.55 X lO^^Re^) 



(10) 



Only for Re > 40 does the growth deviate significantly 
from this form. There is an upturn in the maximum 
growth in approaching Rcc. Above RCc, G"^^^ — oo since 
the flow is linearly unstable and G{t) diverges as r — > c». 
Data up to Re = 46 have been obtained via the tran- 
sient growth calculations described in Sec. [ll] The max- 
imum growth at Rcc is obtained differently. At criti- 
cality, the optimal initial condition is the adjoint eigen- 
mode corresponding to the critical eigenvalue. Under 
linear evolution, this initial condition evolves to the di- 
rect eigenmode. Hence the optimal growth at criticality 
is obtained from the eigenmode and its adjoint. Taking 
llwll = 1 and (m*,m) — 1, i.e. biorthonormalized modes, 
then the maximum growth is given by ||m*||. Based on 
these calculations our estimate of the maximum growth 
at RBc, and hence the maximum growth within the sub- 
critical regime, is G^"^'^ = G'"'''^(Rec) w 6800. 



B. Spatiotemporal Evolution 

We now turn to one of our primary focuses, the tran- 
sient evolution of infinitesimal perturbations. For the 
most part we shall be interested in the qualitative char- 
acter of this evolution and how it depends on Re. 
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FIG. 9. Energy of evolving perturbations computed at Re=40 
for three time horizons: r=20, r=65, and r=95. The three 
curves touch the optimal growth envelope (circles) at their 
respective r values. Qualitatively similar evolution is seen 
over a large range of optimal perturbations. 



We start the perturbation field u' from initial condi- 
tions calculated for optimal linear energy growth and 
evolve the field via the linearized Navier-Stokes equa- 
tions, Eq. (Is]). Recall from Sec. Ill] that the initial con- 
dition giving optimal growth over time horizon r is the 
dominant eigenfunction Vi of .4* (r)^(r). Hence, the ini- 
tial condition and subsequent evolution u'(i), depend on 
time horizon r in defining A* {t)A{t). However, the tran- 
sient dynamics based on a substantial range of r values 
are qualitatively very similar. This is illustrated, in part, 
by Fig. ]9] where we show the the optimal growth envelope 
(denoted by circles) at Re = 40 and also the energy evo- 
lution from optimal initial conditions corresponding to 
three quite different values of r. While there are quanti- 
tative differences between the transient-response curves, 
qualitatively they are similar, suggesting similar fiow fea- 
tures and dynamics are excited by initial conditions op- 
timized across a large range of r values. This holds for 
other values of Re, with the peak in the response curves 
shifting to smaller times for smaller values of Re. In the 
spatiotemporal results which follow, we have opted to fix 
the T at which the optimal perturbations are computed, 
rather than having it vary with Re. All optimal pertur- 
bations are for r=20, as this provides a good choice over 
the whole of the subcritical regime. 

Visualizations of the linear time evolution of perturba- 
tions u'(i), at Re=20 and Re=40 are shown in Figs. 10 



and 
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Perturbation energy, |u'| /2 is plotted in Fig. 



with a fixed energy scale throughout the figure. Figure 
shows the same evolution, except in terms of vorticity. 
Here we visualize not the vorticity in the perturbation 
field itself, but in a superposition of the base flow and 
the perturbation, i.e. u = U -I- eu', where e is chosen so 
that the resulting superposition best resembles what one 
might find in an actual fiow, as for example, might be 



observed in experiments. In all cases only a portion of 
the full computational domain is shown. 

The bottom-most plots correspond to the optimal ini- 
tial conditions u'(i = 0) = Vi. One can see in the energy 
plot that the initial condition is more localized to the 
cylinder at higher Re. In fact the initial condition be- 
comes quite broad spatially at low Re. In the vorticity 
plot one can see the asymmetry of the combined flow 
introduced by the perturbation. The base flow U is sym- 
metric about the centerline, while the perturbation u'(0) 
is antisymmetric. Note, for the Re and r values consid- 
ered here, there is no weak upstream tail in the initial 
conditions of the type shown in Fig. p]^b) , although weak 
upstream tails are found at Re=40 for larger values of r. 
These tails play no qualitative role in the spatiotemporal 
dynamics. 

The perturbation fields are evolved via the linearized 
Navier-Stokes equations, Eq. (Js]), and visualized every 
10 time units. The first obvious point is that in both 
cases the perturbation fields, or more accurately the su- 
perposition of the perturbation fields and base flow, re- 
semble transitory Benard-von Karman vortex streets. At 
Re = 20, the initial perturbation develops into a packet 
of roughly two wavelengths in streamwise extent and ad- 
verts steadily downstream at a speed slightly less than 
1. The peak energy is reached at i ~ 27 and there- 
after the energy decays quite gradually. At Re=40, the 
leading edge of the packet, and the streamwise location 
of the maximum of the response, advects downstream 
at approximately the same speed as at Re=20. In this 
case however, the evolving perturbation develops a long 
trailing series of sinuous oscillations as the excited near- 
wake region undergoes slowly decaying oscillations. The 
streamwise wavelength of oscillations is smaller at Re=40 
than at Re=20. The peak energy at Re=40 is not reached 
until t ~ 60, after the last plot shown. It is evident that 
the growth in the integrated energy of the perturbation 
field is due both to an increase in the maximum point- 
wise energy and also to a significant increase in the spa- 
tial extent of the perturbation field. This second factor 
becomes increasingly important as Re approaches Rcc. 

To further highlight the spatiotemporal character of 
the evolving perturbations, and their dependence on Re, 



we show in Fig. 12 space-time diagrams covering a large 
range of both space and time. Within each tile, space 
is horizontal from a; = to x = 125 and time is ver- 
tical from t = to t = 125. Hence unit speed, that 
of the free-stream velocity, corresponds to 45° in these 
plots. Each row in this figure corresponds to a particu- 
lar Re, from Re=20 to Re=50. The left column shows 
the evolution of energy in the perturbation sampled on 
the flow centerline, that is, contours of u'{x,y = 0,i), 
where u'(t = 0) = Vi is the optimal initial condition. 



For example. Fig. 12 'a) shows the same perturbation as 
in Fig. [Tor a). The contour scale varies from row to row 



and is set so that the maximum energy corresponds to 
white and zero energy corresponds to black. 



The center and right columns of Fig. 12 are explained 
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FIG. 10. (color online). Contours of energy showing the linear evolution of perturbations at Re=20 (left) and Re=40 (right). 
The panels are snapshots at 10 time-unit intervals from t = (bottom) to i = 50 (top). 




FIG. 11. (color online). Same evolution as in Fig. 10 Re=20 (left) and Re=40 (right), but viewed in terms of vorticity. The 



vorticity of a linear superposition of the base flow and the perturbation is shown at snapshots separated by 10 time units from 
t = (bottom) up to i = 50 (top). The maximum vorticity is 1.5. 
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as follows. The center column is the evolution of the 
first sub-dominant optimal mode, that is the evolution 
of u'(i), where u'{t = 0) = V2. The sub-dominant mode 
and subsequent evolution are very similar to that of the 
dominant mode. However, careful inspection shows that 
the sub-dominant mode is spatially phase shifted by a 
quarter wavelength with respect to the dominant mode. 
This is seen as a half wavelength shift in Fig. [12] since a 
pair of vortices (one wavelength) generates two peaks in 
the centerline energy. The pairing of perturbations has 
been observed and discussed elsewhere [231 123 as a man- 
ifestation of streamwise symmetry breaking such that 
modes come in near pairs with similar, but not identical, 
dynamics. The importance of this second, phase-shifted 
mode is that from the pair of modes we can easily con- 
struct an approximate energy envelope eliminating the 
fast oscillations associated with vortex shedding. This is 
shown in the third column where we plot E = Ei + aE2, 
where Ei and E2 are the energy of the dominant (left 
column) and sub-dominant (middle column) perturba- 
tion fields. We choose a so that the peak energy of the 
sub-dominant mode matches that of the dominant mode. 
As one can see this nearly eliminates the fast oscillations 
throughout the space-time plot of E. 

The dynamics seen at Re=20, Re=30, and Re=40 are 
quite similar. There is an increase in energy (both peak 
energy and spatial extent) followed by a decrease with the 
long-term dynamics being a weak wave packet propagat- 
ing and decaying downstream. The effects of varying Re 
in the regime are those already noted: there is a decrease 
in wavelength and an increase near- wake oscillations with 
increasing Re. 

The behavior at Re=45 is, however, qualitatively dif- 
ferent from that seen at Re=40 and below, even though 
Re=45 is still in the subcritical regime. The perturba- 
tion at long times does not have a maximum at some 
downstream location set by how long the perturbation 
has evolved. Instead the maximum is located at a finite 
streamwise location. This is due to the fact that at long 
times the perturbation must evolve to the least stable 
wake eigenmode and there is a qualitative change in the 
spatial structure of this eigenmode at Re ~ 42 (also noted 
by Giannetti and Luchini who give Re ~ 43). Below 
Re = 42 the leading eigenmode is exponentially growing 
downstream and hence appears localized to the down- 
stream computational boundary. Above Re — 42 the 
leading eigenmode has a maximum at finite streamwise 
position with exponential decay far downstream. The lo- 
cation of the maximum decreases as a function of Re and 
is at about x = 34.6 for Re=45. This phenomenon is well- 
known and understood in other systems, e.g. (42H44J . In 
these systems, the switch from downstream growth to 
downstream decay of an eigenfunction occurs when the 
corresponding eigenvalue crosses the essential spectrum. 
The essential spectrum, in turn, is the continuous eigen- 
value spectrum associated with the far-field part of the 
system. It might be of some interest in the future to 
investigate these issues for the cylinder wake. 



For completeness we also show the evolution at Re — 
50, slightly into unstable regime. Although somewhat 
masked by the fact that the perturbation is growing, the 
spatial structure of the mode at long time is not very 
different from that at Re = 45. The perturbation has 
a maximum at about x ~ 19.4 followed by exponential 
downstream decay, matching that of the leading eigen- 
mode from the stability analysis at Re=50. 



C. 3D Energy Growth 

We consider here briefly the energy growth of 3D per- 
turbations, mainly to show that 3D effects are unim- 
portant. The spanwise wavenumber /3 of perturbations, 
Eq. Q, becomes an additional parameter to vary. We 
shall fix the Reynolds number at Re — 40. Optimal 
growth curves over a range of spanwise wavenumbers, 
at representative values of t, are plotted in Fig. [13] and 
growth contours in the /3-t plane are shown in Fig. |14[ 



The thicker line in Fig. 14 denotes the no-growth con- 
tour and energy growth occurs only to the left of this 
contour. Except for small values of r, the growth of 2D 
perturbations (/3 = 0) greatly dominates the growth of 
3D perturbations. 

For short time horizons, (we estimate r < 8.0), the 
largest possible growth is found at nonzero wavenum- 
bers and the range of active wavenumbers increases con- 
siderably as T approaches zero. This shift to high- 
wavenumber modes at short time horizons occurs in other 
shear flows [261 123 145 , but we are unaware of any expla- 
nation for this phenomenon. This does not seem impor- 
tant in any practical sense because the overall response 
of such modes is very small indeed. We have not inves- 
tigated other values of Re in detail, but the unsurprising 
result is that 2D modes dominate the transient response 
in the subcritical wake. 



V. SUMMARY AND DISCUSSION 

We have studied the subcritical response of the cylin- 
der wake by accurately computing the optimal energy 
growth throughout the subcritical regime. We have 
treated at some length the numerical domain require- 
ments for accurate computations within the subcritical 
region. The results themselves show that energy growth 
first occurs as low as Re ~ 2.2, below the onset of separa- 
tion at Re ~ 6.2. Over most of the subcritical regime the 
maximum energy amplification increases approximately 
exponentially in the square of Re. This super-exponential 
dependence on Re is even faster than the exponential de- 
pendence commonly observed in other separated flows 
[23l [26l Hz] ■ However, the maximum growth in the cylin- 
der wake never reaches the extremely large values since 
the wake becomes linearly unstable at a relatively low Re 
where the maximum energy growth is about 6800. 
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FIG. 12. (color online). Plots showing the space-time evolution of energy in the perturbation u'. Each tile covers x between 
and 125 (horizontal) and t between and 125 (vertical), and shows the energy of perturbations sampled on the centerline 
J/ = 0. Each row corresponds to a different Reynolds number, specifically: Re = 20 (a,b,c). Re = 30 (d,e,f), Re = 40 (g,h,i). 
Re = 45 (j,k,l), and Re = 50 (m,n,o). For each Reynolds number we show the evolution of the dominant mode (left column), 
the first sub-dominant mode (center column), and a combination (right column) revealing the envelope of the perturbation 
as explained in the text. The scale of each row of tiles is normalized by the maximum energy in the right column over the 
space-time domain, with white corresponding to the highest energy and black to zero energy. 



We have considered the structure of the optimal tran- 
sient dynamics. The evolving perturbations are of the 
form of transitory Benard-von Karman vortex streets. 
At lower Re wave packets of only a few wavelengths are 
formed which propagate downstream. As Re increases 
the packets extend in length due to the slow decay of 
oscillations in the near wake. At Re ~ 42 the spatial 



structure of the response at long times switches from ex- 
ponentially increasing downstream to exponentially de- 
caying downstream so that at about Re — 42 the response 
at long times has a maximum at a finite streamwise lo- 
cation. Finally, at Re = Rec — 46.6 the wake becomes 
linearly unstable. 

It is of interest to relate our results to the understand- 
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FIG. 13. Optimal growth as function of spanwise wavenumber 
P at Re=40. /3 = is dominant for long time horizons, but 
higher j3 may provide slightly larger growth at short time 
horizons. 




FIG. 14. Contour plot of optimal growth at Re=40. The 
thicker black line denotes the contour of no growth. 



ing of subcritical dynamics arrived at by local stability 
analysis, e.g. [3 |9j [TTJ |46H49| and references therein. In 
brief, from sectional stability analysis of wake profiles, 
the picture of the subcritical region is as follows. Be- 
low Re ~ 5 the wake is everywhere stable. Above Re ~ 5 
there is a region of convective instability behind the cylin- 
der and at Re ^ 25 a pocket of absolute instability ap- 
pears within the region of convective instability. The size 
of the absolute pocket grows with Re and is thought to 
be responsible for the actual instability occurring at Rec, 
although prediction of the transition point has eluded 
local analysis. 

In reality, there are two qualitative changes within 
the subcritical regime: the onset of transient growth 
at Re=2.2 and the switch from downstream growth to 
downstream decay of transient structures at Re ~ 42, 
associated with a corresponding change in the structure 



of eigenmodes. It seems that the first of these, the on- 
set of transient energy growth, could be connected with 
the first appearance of a local convective instability. A 
local pocket of convective instability would indeed cor- 
respond to transient response in a global setting. More- 
over, the Re values for the two event are reasonably close. 
This would appear to corroborate the picture first pro- 
posed by Cossu and Chomaz [21 in the context of the 
Ginzburg-Landau equation. In this picture, one can un- 
derstand the transient energy growth as arising from per- 
turbations traveling through a local region of instability, 
where they are amplified, followed by advection into the 
stable downstream wake, where they decay. We caution, 
however, that the cylinder wake is highly non-parallel in 
the near wake region and it would probably be a mistake 
to connect the transient response and the parallel-flow 
analysis in too much detail. 

There is nothing in the actual transient response cor- 
responding to the local opening of the absolute pocket at 
Re '-^ 25, but neither is there expected to be [IH]. We 
have clearly shown an uneventful evolution of the tran- 
sient response between Re = 20 and Re = 30, and in fact 
up to Re = 40. We are unaware of any local analysis of 
the cylinder wake that predicts the shift from growth to 
decay of modes at Re ~ 42, and this might be interesting 
to investigate in the future. 

There is another way to view the relationship between 
our study and concepts of convective and absolute in- 
stability. This is also closely related to some past and 
ongoing experimental studies [301 ISI] ■ While convective 
and absolute instability are strictly defined for stream- 
wise homogeneous flows, which the cylinder wake is not, 
the change in the linear response at Rec has the essen- 
tial character of the transition from convective to abso- 
lute instability and it commonly referred to using these 



terms. One sees this in our Fig. 12 where the subcritical 
cases, Figs. [I2]^c),[l2]^f),[l2];i) and [12]; 1) have the char- 
acter of convective instability; initial perturbations lead 
to wave packets that advect downstream such that even 
though a perturbation grows (for some time) it is simulta- 
neously advected quickly downstream. The supercritical 



12 o) has the character of absolute instability where 
perturbations grow at fixed streamwise locations. Le- 
Gal and Croquette present nice streakline visualizations 
of the transient wake, qualitatively similar to what is 
shown in our Fig. |11[ and discuss this as evidence of con- 
vective instability in the cylinder wake prior to the onset 
of sustained oscillations. Marias et al. use particle image 
velocimetry (PIV) to obtain more quantitative measures 
of the subcritical response generated by rotary motion 
of the cylinder. In particular they measure front veloc- 
ities and study how these behave as RCc is approached. 
Marias et al. also extract integrated energy from PIV 
data. Transient amplification is indeed observed, fol- 
lowed by exponential decay. However, due to the fact 
that experimental perturbations are introduced by cylin- 
der rotation, and not from the optimal initial conditions 
studied here, quantitative comparisons are not presently 
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possible, but may be pursued in the future. 

Finally, we conclude with the issue of numerical accu- 
racy. Our study has highlighted the importance of en- 
suring the numerical convergence of the computational 
domain. Transient growth problems in open flows with 
inflow-outflow boundary conditions are particularly sus- 
ceptible to deficiencies in the extent of the computational 
domain. This is true not only in the downstream region, 
but in the cross-stream and especially the inflow dimen- 
sions. It is well known that for external flows enforc- 
ing boundary conditions too close to a body can lead to 
deformation of the underlying basic flow [57H40) . Accu- 
rate resolution of perturbation fields for transient growth 
problems can impose yet more severe requirements. The 
cylinder wake is a prime example of a flow in which the 
requisite domain can be far greater for transient growth 
computations than for other types of calculations. 
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FIG. 15. Convergence of the base flow stagnation point with 
mesh dimensions. Points indicate the computed values. In 
(a) Lc is fixed at 25 and in (b) Li is fixed at 25. 



Appendix 

In this appendix we present convergence results for 
base-flow calculations, stability calculations, and poly- 
nomial order. 

Base flow convergence is assessed through two indica- 
tors: the position, Xs, of the stagnation point marking 
the end of the recirculation region and velocity profiles 
just downstream of the cylinder. Figure [15] summarizes 
the convergence of the stagnation point with mesh di- 
mension. The stagnation point is not present at Re = 5, 
and consequently this case does not appear. Percentage 
errors are relative to the calculation using L^ = 65 and 
ic = 65, respectively. The stagnation point is seen to be 
highly converged, as a function of domain dimensions, for 
L, = Lc = 45. 

We may also compare the values of Xg with values re- 
ported in previous studies P^, H-S, 37, 50-53 . Consistent 
with other studies, we find for Re > 6.2 the stagnation 
point obeys 

2;^~0.5 + 0.067(Rc-6.2), 

with specific converged values: Xg = 1.422 at Re = 20 
and Xs = 2.762 at Re = 40. These agree very will with re- 
cent computational studies by Giannetti and Luchini [13] 
and Ye et al. [53] . 

Examination of streamwise velocity profiles is found to 
provide a more detailed view of base-flow distortion due 
to the finite-size effects. Figure [TH] shows velocity profiles 
at location a; = 3. Only a limited cross-stream range in y 



is shown in the vicinity of where the streamwise velocity 
reaches its maximum, as this is where the effects of do- 
main confinement are most pronounced. Constriction of 
the cross-stream mesh leads to an especially inaccurate 
calculation of the base flow, particularly at low Re, while 
the effect of restricted inflow length is less signiflcant in 
general [Fig. 16 a) and (c)]. In any case, the base flow 



is again seen to be highly converged, as a function of 
domain dimensions, for L^ = L^, = 45. 

The dependence of the linear stability calculations on 
domain size is examined through determination of the 
critical Reynolds number, Rec, on different domains. For 
each domain, we compute the base flow and the eigenval- 
ues at Re = 43, 44, 45, and 46. From these we extrapolate 
to find RCc where the real part of the leading eigenvalue 
crosses zero. The results are shown in Fig. |17[ where 
as before we report percentage error in the value of Re^ 
with respect to the value obtained using Li = 65 and 
Lc = 65, respectively. Interestingly, one sees very little 
effect of cross-stream restriction here. In any case, Rcc 
is well determined for Li = L^ = 45, with an error of less 
than 0.1%. 

We may also compare directly the value we obtain for 
Rcc with that obtained in other studies. To three signif- 
icant figures, with Li = Lc = 45, we find 

Rcc = 46.6. 

This value agrees to within half a percent with recent 
stability calculations by Giannetti and Luchini |13) and 
Marquet et al. [25] who quote values of RCc = 46.7 and 
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FIG. 16. Streamwise velocity profiles of base flows at Re=5, 
(a) and (6), and at Re=46, (c) and (d), showing variation 
with Li and Lc- For variations in Li, we fix Lc — 25. For 
variations in Lc, we fix Li = 25. 



Rbc = 46.8, respectively. 

Finally, having set the overall mesh dimensions, we 
consider the convergence of computations with respect to 
the polynomial order of the spectral-element expansion. 
The polynomial order is chosen to ensure there is the 
necessary refinement to resolve the finest characteristics 
of the flow at the highest Reynolds number under con- 
sideration. Base flow and subsequent transient growth 
calculations at Re = 46 have been carried out for a range 
of polynomial orders as summarized in Table |ll] A poly- 
nomial order of P = 6 is found to be sufficient and is 
used for all results reported in Sec. |TV| 
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FIG. 17. Convergence of critical Reynolds number Rcc with 
mesh size. Points indicate the computed values. 



Order 



G{t = 20) 



138.91 
108.29 
156.29 
156.20 
156.19 
156.19 



TABLE II. Convergence of optimal growth results with poly- 
nomial order on the mesh Li — 45, Lc = 45, Lo — 125. The 
base flow and optimal growth G{t = 20) at Re = 46 are both 
computed for the polynomial orders indicated. 
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